home <- here::here()library(tidyverse)library(tidylog)library(RCurl)library(sdmTMB)library(RColorBrewer)library(devtools)library(patchwork)library(ggstats)library(ggh4x)library(viridis)library(sdmTMBextra)# Source map-plot#source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")source(paste0(home, "/R/functions/map-plot.R"))
Reading layer `StatRec_map_Areas_Full_20170124' from data source
`/Users/maxlindmark/Dropbox/Max work/R/cod-interactions/data/shapefiles/ICES-StatRec-mapto-ICES-Areas/StatRec_map_Areas_Full_20170124.shp'
using driver `ESRI Shapefile'
Simple feature collection with 11074 features and 17 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -44 ymin: 36 xmax: 69 ymax: 85
Geodetic CRS: WGS 84
Read data & prepare data
d <-read_csv(paste0(home, "/data/clean/aggregated_stomach_data.csv")) |>drop_na(group) |>drop_na(oxy)t <- d |>filter(fle_kg_km2 >0)t <- d |>filter(mcod_kg_km2 >0)t <- d |>filter(scod_kg_km2 >0)# Calculate relative prey weights (saduria and benthos)d <- d |>rename(oxygen = oxy) %>%mutate(tot_weight =rowSums(select(., ends_with('_tot'))), benthic_weight = amphipoda_tot + bivalvia_tot + gadus_morhua_tot + gobiidae_tot + mysidae_tot + non_bio_tot + other_crustacea_tot + other_tot + other_pisces_tot + platichthys_flesus_tot + polychaeta_tot + saduria_entomon_tot) |>rename(saduria_weight = saduria_entomon_tot,flounder_density = fle_kg_km2,large_cod_density = mcod_kg_km2,small_cod_density = scod_kg_km2) |>mutate(tot_rel_weight = tot_weight / (pred_weight_g - tot_weight), benthic_rel_weight = benthic_weight / (pred_weight_g - tot_weight),saduria_rel_weight = saduria_weight / (pred_weight_g - tot_weight)) |> dplyr::select(-ends_with("_tot")) |> dplyr::select(-predator_latin_name, date) |># Add small constant to large cod density because we want to take the log of itmutate(large_cod_density =ifelse(large_cod_density ==0,min(filter(d, mcod_kg_km2 >0)$mcod_kg_km2)*0.5, large_cod_density),flounder_density =ifelse(flounder_density ==0,min(filter(d, fle_kg_km2 >0)$fle_kg_km2)*0.5, flounder_density),small_cod_density =ifelse(small_cod_density ==0,min(filter(d, scod_kg_km2 >0)$scod_kg_km2)*0.5, small_cod_density)) |># Scale variablesmutate(fyear =as.factor(year),fquarter =as.factor(quarter),fhaul_id =as.factor(haul_id),depth_sc =as.numeric(scale(depth)),oxygen_sc =as.numeric(scale(oxygen)),density_saduria_sc =as.numeric(scale(density_saduria)),flounder_density_sc =as.numeric(scale(log(flounder_density))),large_cod_density_sc =as.numeric(scale(log(large_cod_density))),small_cod_density_sc =as.numeric(scale(log(small_cod_density)))) |># Scale length by group ..mutate(pred_length_cm_sc =as.numeric(scale(pred_length_cm)),.by = group)d |>rename("Relative benthic weight"="benthic_rel_weight","Relative Saduria weight"="saduria_rel_weight") |>pivot_longer(c("Relative benthic weight", "Relative Saduria weight")) |>mutate(group =str_to_sentence(group)) |>ggplot(aes(value)) + ggh4x::facet_grid2(factor(group, levels =c("Flounder", "Small cod", "Large cod"))~name,scales ="free", independent ="all") +geom_density(color =NA, alpha =0.8, fill ="grey30") +labs(y ="Density", x ="Value") +NULL
ggsave(paste0(home, "/figures/supp/data_distribution.pdf"), width =17, height =17, units ="cm")# Size distribution by yeart <- d |>filter(group =="large cod")hist(t$pred_length_cm)
d |>distinct(pred_id, .keep_all =TRUE) |>mutate(group =str_to_sentence(group)) |>ggplot(aes(pred_length_cm, color =as.factor(year))) +facet_wrap(~group, scales ="free", ncol =1) +geom_density(alpha =0.8, fill =NA) +labs(y ="Density", x ="Predator length (cm)", color ="Year") +scale_color_viridis(discrete =TRUE, option ="mako") +NULL
ggsave(paste0(home, "/figures/supp/predator_sizes.pdf"), width =11, height =17, units ="cm")# Sample size per hauld |>summarise(n =length(unique(pred_id)), .by = haul_id) |>ggplot(aes(n)) +geom_histogram()
# How many with 3 or fewer and what are their sample sizes?d |>summarise(n =length(unique(pred_id)), .by = haul_id) |>mutate(s =ifelse(n <3, "s", "l")) |>summarise(n =n(), .by = s)
# A tibble: 2 × 2
s n
<chr> <int>
1 l 316
2 s 29
# (29/316)*100# What is the size range of those 10%?d |>mutate(n =length(unique(pred_id)), .by = haul_id) |>filter(n >1& n <11) |>summarise(max =max(pred_length_cm),min =min(pred_length_cm),.by =c(group, haul_id)) |>mutate(diff = max - min) |>as.data.frame()
group haul_id max min diff
1 large cod 2015_4_6 49 46 3
2 large cod 2015_4_12 49 31 18
3 large cod 2015_4_32 38 27 11
4 small cod 2015_4_32 18 16 2
5 small cod 2015_4_51 21 21 0
6 large cod 2015_4_51 45 45 0
7 large cod 2015_4_53 45 41 4
8 flounder 2015_4_6 39 35 4
9 flounder 2015_4_12 26 22 4
10 flounder 2015_4_32 31 9 22
11 small cod 2016_1_16 18 16 2
12 large cod 2016_1_23 46 29 17
13 large cod 2016_1_37 49 30 19
14 small cod 2016_1_37 25 25 0
15 large cod 2016_1_53 34 31 3
16 large cod 2016_1_73 48 26 22
17 small cod 2016_1_73 22 18 4
18 large cod 2016_1_81 47 47 0
19 small cod 2016_1_81 14 14 0
20 small cod 2016_1_87 17 13 4
21 large cod 2016_1_89 50 48 2
22 small cod 2016_4_17 25 25 0
23 small cod 2016_4_62 14 12 2
24 large cod 2016_4_62 47 47 0
25 flounder 2016_4_17 29 29 0
26 flounder 2016_4_62 35 35 0
27 small cod 2017_1_17 12 11 1
28 large cod 2017_1_60 46 45 1
29 large cod 2017_1_95 45 45 0
30 small cod 2017_1_95 13 12 1
31 flounder 2017_1_95 36 15 21
32 flounder 2017_1_99 38 35 3
33 flounder 2017_1_40 34 20 14
34 small cod 2017_4_52 21 21 0
35 large cod 2017_4_52 28 28 0
36 flounder 2017_4_52 31 20 11
37 large cod 2018_1_52 44 26 18
38 small cod 2018_1_52 25 24 1
39 large cod 2019_4_91 36 29 7
40 flounder 2019_4_82 24 21 3
41 flounder 2019_4_91 26 21 5
42 large cod 2020_1_65 27 27 0
43 flounder 2020_1_57 23 15 8
44 flounder 2020_1_65 26 21 5
45 small cod 2020_4_213 24 24 0
46 flounder 2020_4_213 28 20 8
47 flounder 2020_4_214 27 22 5
48 large cod 2021_1_56 33 26 7
49 small cod 2021_1_56 23 20 3
50 large cod 2021_1_83 28 26 2
51 large cod 2021_1_84 31 31 0
52 small cod 2021_1_84 22 10 12
53 small cod 2021_1_85 12 10 2
54 large cod 2021_1_86 32 29 3
55 small cod 2021_1_86 10 10 0
56 large cod 2021_1_89 29 27 2
57 large cod 2021_1_91 37 37 0
58 small cod 2021_1_91 24 21 3
59 large cod 2021_1_98 40 27 13
60 small cod 2021_1_98 23 23 0
61 large cod 2021_1_99 40 27 13
62 small cod 2021_1_99 25 19 6
63 large cod 2021_1_100 38 38 0
64 small cod 2021_1_100 20 20 0
65 large cod 2021_1_103 44 38 6
66 small cod 2021_1_103 25 10 15
67 large cod 2021_1_107 36 33 3
68 large cod 2021_1_108 41 31 10
69 small cod 2021_1_108 10 10 0
70 large cod 2021_4_258 34 34 0
71 small cod 2021_4_258 23 23 0
72 flounder 2021_4_258 25 18 7
73 flounder 2022_1_88 26 22 4
74 large cod 2022_4_241 32 32 0
75 large cod 2022_4_266 40 29 11
76 small cod 2022_4_266 23 19 4
77 large cod 2022_4_275 34 29 5
78 flounder 2022_4_241 34 20 14
79 flounder 2022_4_266 24 23 1
#summarise(mean = mean(diff))
Quick explore
Correlation between variables
# Plot correlation between variablesd_cor <- d |> dplyr::select("oxygen_sc", "density_saduria_sc", "flounder_density_sc","large_cod_density_sc", "small_cod_density_sc", "depth_sc")corr <-round(cor(d_cor), 1)library(ggcorrplot)ggcorrplot(corr, hc.order =TRUE, type ="lower", lab =TRUE, lab_size =2.5) +theme_classic() +theme(axis.text.x =element_text(angle =90, vjust =0.3))
Sample size
d |>group_by(species) |>summarise(n =n())
group_by: one grouping variable (species)
summarise: now 2 rows and 2 columns, ungrouped
# A tibble: 2 × 2
species n
<chr> <int>
1 Cod 5484
2 Flounder 3851
d |>group_by(species, quarter) |>summarise(n =n())
group_by: 2 grouping variables (species, quarter)
summarise: now 4 rows and 3 columns, one group variable remaining (species)
# A tibble: 4 × 3
# Groups: species [2]
species quarter n
<chr> <dbl> <int>
1 Cod 1 2906
2 Cod 4 2578
3 Flounder 1 2081
4 Flounder 4 1770
Fit models
Groups are: small cod, large cod and flounder. Response variables are: saduria_rel_weight, benthic_rel_weight or total weight. The latter is only for adult cod, because essentially all prey are benthic for small cod and flounder.
Model random effect structure is selected with AIC (see script 02-)
# This is the reason we don't do total weight for flounder and small codd |>filter(tot_rel_weight >0) |>group_by(group) |>mutate(ben_prop = benthic_rel_weight / tot_rel_weight) |>summarise(mean_ben_prop =mean(ben_prop))
# A tibble: 3 × 2
group mean_ben_prop
<chr> <dbl>
1 flounder 0.978
2 large cod 0.591
3 small cod 0.956
Covariates are: ~ 0 + fyear + fquarter + depth_sc + spatial + spatiotemporal random fields + density covariates. For saduria, we use saduria also in interaction with cod and flounder. For cod we use small cod because large and small cod are very correlated. For benthic and total prey, we instead use oxygen, more as a proxy, as the interaction variable
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
#pal <- brewer.pal(n = 8, name = "Dark2")[c(2, 7, 6)]pal <- (brewer.pal(n =11, name ="RdYlBu")[c(11, 4, 1)])range_df |>arrange(estimate)
# A tibble: 7 × 7
term estimate std.error conf.low conf.high group model
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 range 9.59 2.52 5.74 16.1 large cod benthos
2 range 18.5 4.39 11.6 29.5 flounder benthos
3 range 19.5 5.59 11.2 34.2 small cod benthos
4 range 25.6 6.55 15.5 42.3 large cod total
5 range 29.8 10.9 14.5 61.2 large cod saduria
6 range 31.7 10.3 16.8 59.8 small cod saduria
7 range 72.5 16.4 46.5 113. flounder saduria